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Inverse methods of statistical mechanics are becoming productive tools in the design of materials 
with specific microstructures or properties. While initial studies have focused on solid-state design 
targets (e.g, assembly of colloidal superlattices), one can alternatively design fluid states with desired 
morphologies. This work addresses the latter and demonstrates how a simple iterative Boltzmann 
inversion strategy can be used to determine the isotropic pair potential that reproduces the radial 
distribution function of a fluid of amorphous clusters with prescribed size. The inverse designed pair 
potential of this “ideal” cluster fluid, with its broad attractive well and narrow repulsive barrier at 
larger separations, is qualitatively different from the so-called SALR form most commonly associated 
with equilibrium cluster formation in colloids, which features short-range attractive (SA) and long- 
range repulsive (LR) contributions. These differences reflect alternative mechanisms for promoting 
cluster formation with an isotropic pair potential, and they in turn produce structured fluids with 
qualitatively different static and dynamic properties. Specifically, equilibrium simulations show that 
the amorphous clusters resulting from the inverse designed potentials display more uniformity in 
size and shape, and they also show greater spatial and temporal resolution than those resulting from 
SALR interactions. 


1. INTRODUCTION 

The computational design of interactions for targeted 
self assembly is a powerful approach in the search for 
new materials with specified microstructures, properties, 
or functionality. It is typically pursued via a strategy 
where the macroscopic behaviors of a subset of promis¬ 
ing systems with different microscopic interactions (e.g., 
patchiness [1-6] or shape [6-8]) are characterized by ex¬ 
tensive “forward” molecular simulation calculations and 
compared to one another using appropriate figures of 
merit. Such forward approaches have been instrumental 
in discovering novel organizational motifs in crystalline 
or microphase-separated solids. However, since forward 
calculations (or experiments) can be expensive and time- 
consuming for complex systems, this method is perhaps 
most useful where physical intuition can guide the se¬ 
lection of the microscopic interactions to be considered. 
For example, possible locations and sizes of attractive 
“patches” on colloidal particles could be chosen a priori 
by considering how these variables affect mutual patch 
alignment with nearest neighbors when particles are in 
a desired superlattice structure versus other competing 
morphologies [2, 3, 9-13]. 

For less intuitive materials design problems, system¬ 
atic alternatives to forward searches may be helpful. In¬ 
verse methods of statistical mechanics, which formally 
optimize microscopic interactions toward attainment of 
a desired macroscopic outcome, are one such emerg¬ 
ing class of complementary techniques [14-17]. Inverse 
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approaches have been recently applied to gain insights 
into nontrivial materials design problems including the 
search for isotropic, repulsive interactions that can sta¬ 
bilize low-coordinated crystals in two (e.g., honeycomb 
and square) [18-20] or three (e.g., diamond and simple 
cubic) [19, 21-23] dimensions. However, such methods 
are not limited to designing solid-state targets, and stan¬ 
dard tools could in principle be exploited to find inter¬ 
actions that imbue equilibrium fluid states with desired 
microstructural features. Some of the organizational mo¬ 
tifs of these designer fluids could, in turn, be captured in 
non-equilibrium solid states (e.g., gels or glasses formed 
from the flud via a rapid quench or compression). 

Here, we present a simple methodology for the inverse 
design of fluid structure via optimization of an isotropic 
pair interaction. It comprises two steps: generation of 
a configurational ensemble of target microstructures via 
simulations using an artificially complex, many-body in¬ 
teraction chosen to guarantee assembly of the desired 
morphology, and then use of a tool from systematic 
coarse-graining [24-29] to reduce the many-body inter¬ 
action to an effective pair potential. In the present work, 
we adopt iterative Boltzmann inversion (IBI) for the lat¬ 
ter step, which uniquely determines the pair potential 
that will generate the radial distribution function (RDF) 
of the target ensemble at equilbrium. As a first applica¬ 
tion of this methodology, we attempt to inverse design 
a pairwise potential that forms a fluid of “ideal” amor¬ 
phous equilibrium clusters of prescribed size. Clustered 
fluids of colloidal particles have attracted considerable 
interest due their novel multiscale structure, their rich 
dynamic and rheological properties, and their potential 
functionality [30-45]. 

The classic paradigm for forming equilibrium clusters 
from an isotropic pair potential focuses on models that 
exhibit a combination of short-range attractive (SA) and 
longer-range repulsive (LR) contributions, commonly re- 
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ferred to as an SALR model [30-38, 44, 45]. Various 
functional forms for the attractions have been studied 
(modeling, e.g., polymer-mediated depletion forces be¬ 
tween colloids), typically in combination with a repulsive 
Yukawa tail to model weakly screened Coulombic inter¬ 
actions. The attractions drive particle association, but 
the longer-ranged replusions lead to self-limited growth 
(i.e., finite sized aggregates). In contrast to systems lack¬ 
ing competitive repulsions, the formation of clusters in 
SALR fluids can either suppress macroscopic phase sep¬ 
aration to lower temperature and density or eliminate it 
altogether [30, 31, 46]. 

Although the SALR model qualitatively captures the 
effective pair potential and equilibrium cluster behavior 
seen in some types of experimental systems (e.g., mix¬ 
tures of charged-stabilized colloids with weakly interact¬ 
ing polymers), it does not generate microstructures re¬ 
flecting properties typically expected in the idealized pic¬ 
ture [47] of such cluster phases [31 38, 45, 46] that are 
found in other experimental systems [48, 49], which we 
denote here as ideal cluster (IC) fluids: particle assem¬ 
blies that are monodisperse, spherical, long-lived, and 
fluid-like in terms of inter- and intra-cluster structure and 
mobility. Recently, Glotzer, Kotov, and coworkers have 
demonstated that a many-body potential which incor¬ 
porates environment-dependent charge renormalization 
during assembly can lead to clusters that are monodis¬ 
perse, spherical, and amorphous [48-50]. The role of sur¬ 
face charge renormalization has been studied by others as 
well [51, 52]. However, whether many-body interactions 
are in general neccessary to assemble the complex multi¬ 
scale structures of IC fluids has been an open question. 
Additionally, keeping to the level of a pair interaction 
has the benefit that, through a systematic mapping, it 
can be regarded as a low dimensional approximation of 
a given many-body interaction. The simplified form can 
then yield key physical insights. 

Here, we show that one can inverse design pair poten¬ 
tials that readily assemble into IC fluids under equilib¬ 
rium conditions. Interestingly, these potentials exhibit 
a broad attractive well together with a narrow repulsive 
barrier at larger separations, which-while also a competi¬ 
tive balance between two interactions-is qualitatively dif¬ 
ferent from those of the SALR fluid. These differences 
imply distinct physics governing cluster formation in IC 
and SALR fluids, and we compare the static and dynamic 
properties of clusters in these systems, introducing a new 
metric for cluster lifetime to quantitatively characterize 
the latter. In the analysis, we also discuss practical as¬ 
pects in the inverse design of pair potentials for complex 
fluids using IBI. 

The remainder of this paper is organized as follows. 
Section 2 outlines the constrained model used to gener¬ 
ate configurational ensembles of the targeted ICs, the IBI 
inverse design method employed to discover the final IC 
potentials, the SALR model used and the metric we in¬ 
troduce for the cluster lifetime analysis. Results of the 
inverse cluster design and comparisons beween IC and 


SALR fluids are presented in Section 3. The paper is 
concluded in Section 4 with a discussion of future goals 
and possible improvements to the approach. 


2. METHODS 

2.1. Constrained Monte Carlo simulations 

The first step in the inverse design of ICs is to produce 
a physically realizable target RDF, gt g t{r), correspond¬ 
ing to a configurational ensemble with the desired struc¬ 
tural properties. To generate such RDFs, we utilize con¬ 
strained Monte Carlo simulations of N = 2048 total hard 
core (HC) particles of diameter d , which are divided into 
equisized amorphous assemblies of either nt g t = 8, 16, or 
32 particles, each representing a single cluster. To en¬ 
force cluster association, single-particle translations are 
constrained by a many-body intra-cluster potential act¬ 
ing on the instantaneous cluster radius of gyration, R, 
as 

/3^ intra (R) = A(R 2 -R 2 ) 2 (1) 

where /3 = l/(kBT) (fee is Boltzmann’s constant and T 
is temperature), A is a positive scalar amplitude, and R 
is a target radius of gyration. For a given nt g t there is a 
practical lower limit to what R values can be sampled by 
a cluster due to hard-core packing constraints; thus, any 
R below this limit yields virtually identical behavior for 
appropriately chosen values of A. For ?r tg t = 8, 16, and 
32, we use R = 0.6d, 0.8 d, and 1.2d and A = 300, 265, 
and 170 respectively. For all particle packing fractions 
?7 = (n/6)Nd 3 /V studied (where V is volume), these 
parameters yield corresponding average radii of gyration, 
(R) « 0.860(1, 1.105d, and 1.476d. The insensitivity of 
(R) with respect to ry for a given nt g t is due to the preset 
compactness of the clusters. 

We also introduce a longer-ranged Yukawa repulsion 
between the cluster center-of-mass (COM) pairs, which 
improves convergence of the IBI scheme (discussed be¬ 
low) and is defined by 

/3<pcoM(rcoM) = -exp[—rcoM/z] (2) 

rcoM 

where rqoM is the pair COM distance between two clus¬ 
ters and B and 2 are the repulsive amplitude and range, 
respectively. We set z = 0.12 for all systems and B = 
1 x 10 12 , 1 x 10 15 , and 1 x 10 21 for nt g t = 8, 16, and 32, re¬ 
spectively, for all ?y. These parameters furnish very steep, 
hard-core-like repulsions around the clusters with effec¬ 
tive hard-core diameters of d e ff ~ 3.18d, 3.98 d, and 5.60d, 
respectively. From d e s, we can also obtain the effective 
volume fraction of whole clusters (treating them as renor¬ 
malized objects) via the expression ?y e ff = rjd 3 s /{n tg td 3 ). 

Given these definitions, we propagate the Monte Carlo 
trajectories via cluster COM translational moves (10%) 
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and single-particle displacements (90%). Note that clus¬ 
ter rotational moves are unnecessary as the single-particle 
moves are sufficient to randomize the intra-cluster struc¬ 
tures. 

Once gtgt(r) has been obtained, we must smooth out 
the discontinuous peak at contact that results from the 
use of the hard-core constraint, so as to be consistent 
with the use of continuous pair potentials as required 
by the molecular dynamics IBI framework. To do this, 
we construct an approximate hard-core mapping onto a 
steep, purely repulsive 50-25 Weeks-Chandler-Andersen 
(WCA) potential [53] (in dimensionless form) 


PvwcA(r) = H(r 0 -r) 


d-WCA 


50 


dwcA 


25-1 


+1 


( 3 ) 

where H(x) is the Heaviside step function, dwCA is the 
effective core diameter and ro = 2 1 / 25 dwcA)- The map¬ 
ping uses a linear extrapolation of the hardcore target 
RDF, <7tgt,Hc(j")> near contact into the core region to 
locally approximate the well known cavity distribution 
function, j/tgt,Hc(0 [53]. It is generally accepted that 
ytgt.HcW « ytgt.scM, where 2/ tg t,sc(0 is any steeply 
repulsive soft-core (SC) analog that approximates the 
hard core version via some non-unique mapping crite¬ 
rion. For our purposes, the simplest mapping, dwCA = d, 
is sufficient. The final soft-core profile is constructed as 
9t g t,wCA(r) ss exp[-/3<p W cA(?’)]ytgt,Hc(0, which we de¬ 
note as gtgt{r) in the remaining sections. 


2.2. Iterative Boltzmann inversion 


IBI is a conceptually simple and popular approach for 
solving the inverse statistical-mechanical problem of dis¬ 
covering the unique pair-potential u(r ) corresponding to 
a particular RDF [24-28]. In general, there is no guar¬ 
antee that such a potential exists according to the Hen¬ 
derson theorem [54]; however, if the potential exists, IBI 
is a suitable tool for recovering it. Inverse designed po¬ 
tentials depend on the state point of interest (?? and T 
dependent); however, varying T at fixed ij leads to trivial 
rescaling of the potential, thus all potentials are reported 
in units of thermal energy for generality. The explicit 
density dependence of our potentials is discussed in Sec¬ 
tion 3.2. 

The IBI procedure requires an initial-guess potential 
Mi(r), which at the lowest densities we take to be the 
target potential of mean force u±(r) = — fcsTlnfgtgtM]- 
At higher densities, we use converged results from the 
lower densities. Simulation of u\(r) provides the first 
trial RDF, gi(r), and a new potential is calculated ac¬ 
cording to the general formula 


it i+ i(r) = Ui{r) + a m k B T In 


' 9i(r) ' 
.fftst(r). 


( 4 ) 


where a m is a mixing parameter to help control the con¬ 
vergence. The simulation step and potential update steps 


are carried out successively until satisfactory convergence 
in u(?’) [g(r)] is achieved. For our highly structured RDFs 
(orders of magnitude variation; see Fig. 1), a m is best 
kept very small and a m ~ 0.005 — 0.02 provides the best 
convergence while maintaining a quasi-equilibrium sys¬ 
tem during the iterative procedure. In practice, the po¬ 
tential is also always cut and shifted at a lengthscale 
r c after each IBI iteration. For the work here, a value 
of r c = 8 d was sufficient and could be lowered in some 
cases, e.g., when considering smaller clusters. 

Implementation of IBI is accomplished through the 
Versatile Object-oriented Toolkit for Coarse-graining Ap¬ 
plications (VOTCA) [28], which is implemented with the 
GROMACS 4.5.3 molecular dynamics (MD) package [55]. 
We perform simulations comprising N = 2048 pa rticles 
using a time step of dt ss 0.0003y / d 2 r?z/(fc_BT) (m is 
the particle mass) and a velocity-rescale thermostat for 
T with characteristic time constant r = lOOdt, where 
rescaling is done every 10 dt. VOTCA utilizes GRO¬ 
MACS trajectories to calculate RDFs and potential up¬ 
dates accoring to Equation 4. 


2.3. SALR model systems 

To contextualize the behaviors of the newly designed 
IC potentials, it is useful to make comparisons to results 
of an SALR interaction potential known to exhibit equi¬ 
librium cluster phases [30-38]. Specifically, we compare 
IC results with those from a ternary mixture model de¬ 
veloped in a previous study [31] that can generate both 
amorphous and microcrystalline clusters. The pair po¬ 
tential in this model is defined as 

_* I £ 

P<PsL\i,j{xi,j) = 4[x+(l-2^, J )A x ](a;r 2a -xr j a )+Q e 

x i,j / S 

.( 5 ) 

where x = r/d is a non-dimensionalized interparticle 
separation, d is the characteristic particle diameter, x 
quantifies a short-range attractive strength (we choose 
a = 100 to set attractive wells of 0(1%) of the core 
diameter), and Q and £ respectively set the magnitude 
and range of a long-range Yukawa repulsion. Addition¬ 
ally, Sij is the Kronecker delta, with i (or j) = —1,0,1 
corresponding to small, medium, and large particles, re¬ 
spectively; the generalized interparticle distance is de¬ 
fined Xij = x — (l/2)(* + j)(Ad/d). The remaining pa¬ 
rameters are perturbative shifts to particle size Ad and 
energy A x , respectively. The values of x, Q , and £ were 
tuned to generate SALR clusters of comparable size to 
the optimized ICs. 

To generate amorphous cluster phases, we follow our 
previous work [31] and examine three-component mix¬ 
tures that approximate suspensions with 10% size poly- 
dispersity by using 20% small, 60% medium, and 20% 
large particles with size perturbation = 0.158d, where 
we set A x = 0.25x to gently promote mixing. To ex¬ 
amine the more commonly studied monodisperse single- 
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FIG. 1. Comparison of the target (hard-core) and optimized 
radial distribution functions g(r) for ntgt = 8, 16, and 32 at 
packing fraction r/ = 0.04. 


on the presence of a local maxima in the CSD at n* oc¬ 
curring in the range 1 -C n* -C N. 

Calculation of CSDs, which depend on many-body in¬ 
terparticle correlations, provides a practical means to ac¬ 
cess important information regarding how well clustering 
is reproduced in our pair-potential system. In particular, 
matching RDFs via IBI does not guarantee that the cor¬ 
rect cluster size is reproduced, nor does it provide infor¬ 
mation on any undesired (i) polydispersity of emergent 
aggregates and (ii) free monomer or other small aggre¬ 
gates, both of which are absent in our constrained, target 
simulations. 

As a complement to the CSD analysis, we use a pre¬ 
viously published approach [57] to calculate probability 
distributions P(x) where x = < 74 , qe , W 4 , and wq are 
the four standard parameters characterizing local bond- 
orientational order (BO). [58]. The comparison of these 
distributions between the constrained and optimized sys¬ 
tems provides a first-order quantification of how effec¬ 
tively the RDF mapping preserves higher-order, local 
structural correlations within clusters. In calculating the 
local BO parameters, we employ the same r cu t used in 
the CSD calculation for identification of each particle’s 
nearest neighbors. 


component model that exhibits microcrystalline clusters, 
we simply set = A x = 0. 

With the SALR model, we perform three-dimensional 
MD simulations of N = 2960 particles in the canoni¬ 
cal ensemble with periodic boundary conditions using 
LAMMPS [56]. We use an integration time-step of 
dt = 0.0005 ^d 2 m/(kBT) 1 include interactions out to a 
cut-off distance of r c = 8 d, and fix temperature via a 
Nose-Hoover thermostat with time-constant r = 2000dt. 


2.4. Cluster-size distribution and bond order 
analysis 

To characterize the instantaneous scale of the equi¬ 
librium cluster aggregates, we calculate cluster-size dis¬ 
tributions (CSDs) quantifying the probability P(n) of 
observing aggregates comprising n particles. As is cus¬ 
tomary [30-38], two particles are considered part of the 
same cluster if at least one of the following conditions 
are met: ( 1 ) their centers are within a pre-defined cut¬ 
off distance r cu t, making them direct neighbors; and/or 
( 2 ) their centers are both within r cut of one (or more) 
other particle(s), i.e., are connected via some perco¬ 
lating pathway. For the IC systems, we generally use 
r cu t = 1.25d (as discussed later, results are not very sen¬ 
sitive to the choice of r cu t); for the SALR systems, we 
choose r cu t to be the range of the attractive well, which 
varies slightly depending on choices of parameters but is 
generally around r cut — 1.05d. Throughout the remain¬ 
der of the manuscript, a phase is considered clustered 
with aggregates of a preferred equilibrium size n* based 


2.5. Cluster persistence 


Despite the frequent measurement of static CSDs, 
there is little discussion in the literature regarding the 
dynamic stability (“lifetime”) of the contributing aggre¬ 
gates. In fact, we are not aware of any generally accepted 
methods that quantify cluster dynamics in particle-based 
systems other than measurements of monomer mean- 
squared displacements [59] and dynamic structure fac¬ 
tors [60], especially in ways that captures cluster persis¬ 
tence by explicitly incorporating dynamic bond informa¬ 
tion. To facilitate such measurements, we introduce the 
correlation function <f>(f), which quantifies the fractional 
similarity of associates in clusters (FSAC) at an initial 
time t = 0 and a lag-time t > 0 , where particles are as¬ 
sociates at a given time if they are in the same cluster. 
More specifically, the correlation function can be written 


N 




i= 1 


shared if) 
total (t) 


( 6 ) 


where db, shared!/) counts the the number of particle i 
associates that are common to t = 0 and t > 0 while 
‘f’i,total (t) is the combined sum of particle i associates 
at t = 0 and t > 0 (without double counting particles 
common to both times). 

For example, if particle 1 is in a cluster with particles 
{2, 3, 4} at t = 0, and in a cluster with particles {2, 3, 5} 
at some t > 0 , then 'bqtotaKi) = 4 while $1 .shared if) 2 . 
The fractional similarity of associates to particle 1 is then 
$i,shared(0/ $ i,total(i) = 0.5. In the special case that 
particle i is a monomer at both time-points we assume 
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$i,shared (£)/$*,total (i) = 1 , thus ensuring monomers that 
remain as monomers contribute positively to the score. 
Averaging the score of every particle then yields Equation 
6 . 

The correlation function has the range $(<) = [0,1], 
where $(t > 0 ) = 1 means that all clusters (includ¬ 
ing monomers) contain the same particles at both time- 
points and = 0 means that all particles possess tem¬ 
porally exclusive sets of associates. In a macroscopic 
system with finite-sized clusters, <&(f —» oo) = 0 due 
to single-particle and cluster diffusion. However, this is 
not captured in a finite-sized box as particles can wrap 
through the periodic boundaries. As such, we include an 
additional rule: if two particles i and j ever move further 
than a half box-length away from one another at t > 0 , 
each is subsequently treated as a new, previously unrec¬ 
ognized particle with respect to the other at all future 
time-points t' > t. Thus, two associates at t = 0 that 
diffuse very far from one another before becoming asso¬ 
ciates again (in a potentially arbitrary cluster) are not 
counted as temporally common. 

Altogether, the counts in Equation 6 can be expressed 
more formally as ^ shared if) — 

and $i,total(*) = X) j ^ i [0ii(O)+0 i j(t)]-$i, s hared(i) where 


By (t) 


1 , i and j in same cluster at lag-time t 
0 , otherwise 


where the instantaneous cluster analysis at each time is 
done in the usual way. The factor 7 ^ (t) enforces the rule 
concerning temporal pairwise drift and is given by 


Vt ' < < i/2 

| 0 , otherwise 

where L is the simulation box length and Vij(t) = 
||r^-(0)-|-Ar^(t)|| with r™ (0) the wrapped initial displace¬ 
ment vector between particles i and j and Ar^- (t) the cor¬ 
responding net cumulative unwrapped displacement over 
lag-time t. 


3. RESULTS AND DISCUSSION 
3.1. Cluster structure 

In Fig. 1, we compare RDFs obtained from the con¬ 
strained Monte Carlo simulations to those of MD sim¬ 
ulations using the IBI-optimized potentials. The g(r) 
profiles show very good agreement over three orders of 
magnitude (demonstrating the successful application of 
the IBI approach), and also-as expected-exhibit features 
that are consistent with clustering and atypical of simple 
fluids. One such feature is the highly structured, liquid 
droplet envelope extending over multiple particle diam¬ 
eters. This droplet region is terminated by a particle- 


rarefied window, which, in a highly averaged sense, de¬ 
fines the cluster center-to-surface distance and provides 
an intuitive division of <? tg t(r) into intracluster and in¬ 
tercluster particle correlations. Intercluster gt s t(i~) corre¬ 
lations are oscillatory, with relatively long characteristic 
wavelengths set by the effective cluster size d e g, which is 
originally encoded during the reference constained MC 
simulations (see Section 2.1). Importantly, d e g also con¬ 
trols the depth of the depletion window in £/tgt(f), where 
the depletion depth is greater when the repulsive-shell 
lengthscale (d e g) is larger. Intercluster repulsion at the 
COM level is essential towards achieving convergence in 
the IBI scheme, as in its absence, we find that the inter¬ 
mediate IBI steps become unstable towards large scale 
aggregation. In addition, the thickness of the repulsive 
layer needed for convergence grows with cluster size (see 
Section 3.4). 


V 

0.02 

0.04 

0.06 

0.08 

n* 

8 

8 

8 


(n) 

7.95 

8.26 

8.48 


5n 

2 

2 

2 


n* 

16 

16 

16 

16 

(n) 

15.90 

16.03 

16.65 

16.25 

8n 

3 

3 

3 

3 

n* 

33 

32 

33 


(n) 

32.68 

31.46 

32.66 


5n 

4 

4 

5 



TABLE I. Average cluster size and polydispersity measures 
for all cluster sizes (ntgt = 8, 16 and 32) and volume fractions. 

Of course, convergence in g(r) alone does not guarantee 
preservation of higher-order particle correlations, which 
could [61] play an important role in clustering behavior. 
However, we find that strong clustering emerges using the 
designed pairwise potentials, as evidenced by the repre¬ 
sentative simulation snapshots shown in Figs. 2(a-c), the 
three corresponding CSDs in Figs. 2(d-f), and simulation 
movies of these systems (provided in the Supplementary 
Material). Visual inspection of Figs. 2(a-c) reveals that 
the self-assembled clusters are ( 1 ) highly-distinguishable 
(well-defined); (2) spherical and droplet-like; and (3) sim¬ 
ilar in size to the enforced analogs in Fig. 1. The CSDs, 
calculated using r cut = 1.25d[62], confirm that the opti¬ 
mized potentials promote clusters of the desired sizes, as 
quantified by both the characteristic maximum n* 

P{n*) = max„P(n) (7) 

and the average value 

OO 

(n) = nP(n) ( 8 ) 

71=1 

both listed in Table 1. We also note that the clusters are 
so well defined that peaks corresponding to infrequently 
connected clusters (two times primary cluster size) are 
also observed at higher 77 (not visible for n tg t =32). 
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FIG. 2. (a-c) Simulation snapshots for ntgt = 8 , 16, and 32, respectively, at a volume fraction of 77 = 0.06. (d-f) Cluster size 
distributions for ntgt = 8 , 16, and 32, respectively, for various 77 . 


To complement the size measures n* and (n), we also 
calculate the peak-width Sn according to 

«tgt+<Sn 

0.90 < p ( n ) ( 9 ) 

n=ntgt— Sn 

which is the (one-sided) range in n about n* that ac¬ 
counts for 90% of the summated P{n) curve. As demon¬ 
strated in Table 1, Sn is of a reasonable size for all ?r tgt , 
with the ratio Sn/n* decreasing with increasing n* (in¬ 
tuitively, Sn/n* —> 0 in the thermodynamic, cluster-size 
limit). Finally, as is clear from Figs. 2(d-f), these highly 
monodisperse clusters also coexist with a small but nu¬ 
merically detectable fraction of free monomer, which can 
also be visually gleaned from Figs. 2(a-c) and the three 
corresponding supplemental movies. This bifurcation 
of the system into two primary populations suggests a 
similarity to liquid-vapor coexistence, albeit on the mi¬ 
croscale. 

While the CSD calculations clearly demonstrate that 
target cluster sizes are reproduced and that minimal in¬ 
trusive smaller or larger objects are present, we also 


find that reproducing pair structure carries over to some 
higher-order structural correlations. Here, we specifi¬ 
cally consider the local higher-order bond-orientational 
order (BO) parameters [57, 58] < 74 , q 6 , uq, and wq, which 
are understood to be strongly correlated with pair struc¬ 
ture for homogeneous disordered liquids and even glassy 
packings [63, 64]. In Fig. 3, we compare probability dis¬ 
tributions P{x) of these local BO parameters for con¬ 
strained and optimized systems of ntgt = 32 clusters at 
77 = 0.06. For the latter case, a distribution of cluster 
sizes is present. Therefore, to compare clusters of simi¬ 
lar size, we collect statistics over clusters of instantaneous 
size n*ztSn (see Table 1) from the optimized simulations. 
The intracluster higher-order correlations are indeed well 
preserved between two cases, though we note that the ob¬ 
served agreement slightly diminishes as target cluster size 
decreases (see Supplementary Material for n tg t = 8 and 
16 calculations), a trend likely related to the issues we 
discuss in the following paragraph. 

Beyond the cluster sizes (8 < n tg t < 32) considered 
in Figs. 1- 2, we also attempted the IBI approach to ob¬ 
tain pair potentials it( 7 ') that would generate smaller and 
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FIG. 3. Bond-orientational order probability distributions for 
ntgt = 32 at packing fraction 77 = 0.06. The filled grey and 
open red curves indicate constrained and optimized simula¬ 
tions (using clusters with sizes between n* ±<5n), respectively. 
Analogous calculations for the two other cluster sizes are avail¬ 
able in the Supplementary Material. 


larger amorphous clusters; however, various challenges 
emerge in either limit. For clusters of size ntgt < 8 , 
we find that it is difficult for an isotropic pair potential 
to generate well-differentiated amorphous clusters while 
simultaneously suppressing similarly sized, ordered , intr¬ 
acluster configurations that cannot be easily penalized 
on the basis of single length-scale. This type of issue 
arises in a minor way even for nt g t = 8 , as evidenced by 
the small n = 13 peak in the CSDs of Fig. 2(d), which 
corresponds to clusters with a central seed and 12 closely- 
packed neighbors (the sphere kissing number [65]). We 
found that the intrusion of ordered off-target clusters 
is most prevalent when attempting to stabilize dimers 
(fitgt = 2), where 3-mers and 4-mers were always more 
highly favored. This is a result of the minimalistic fluc¬ 
tuation (single particle) being of order the size of the de¬ 
sired cluster (dimer) and the fact that 3-mers and 4-mers 
can assemble into objects with no discriminatory length 
scales (triangles and pyramids respectively) for a pair po¬ 
tential to disfavor. More generally though, larger (than 
the target size) intrusive clusters must be ordered and 
closely packed so as to utilize space more efficiently and 
avoid sampling the growth limiting repulsions present in 
the potentials (discussed in Section 3.4) optimized for the 
less space efficient amorphous clusters. 

On the other hand, attempts to generate potentials 
that stabilize fluids of n tg t = 64 failed to converge for all 
77 > 0.02. It is unclear whether ICs could be designed 
for clusters of this size with a different choice of param¬ 


eters in the target simulation (e.g., we found that in¬ 
creasing d e ff resulted in a somewhat more stable, though 
ultimately unsuccessful, optimization) or if this failure is 
symptomatic of a fundamental limitation of a pair po¬ 
tential to stabilize fluids of large ICs. Exploring and ar¬ 
ticulating the limits of a pair potential to create given 
fluid architectures remains an open question for future 
research. 


3.2. Dependence on density 

While the results above demonstrate that the IBI op¬ 
timization generates pair potentials that induce the de¬ 
sired clustering, we next demonstrate that these poten¬ 
tials are also robust with respect to changes in 77 via two 
complementary approaches. The first is a comparison 
of the potentials optimized at specific values of 77 in the 
supplementary material. Overall insensitivity, including 
the functional form, to density is found across all clus¬ 
ter sizes-only a weak decrease in the overall amplitude 
of the potentials with increasing 77 is observed. The sec¬ 
ond demonstration of insensitivity to density is demon¬ 
strated by simulating optimized potentials generated for 
a particular 77 under either a slow (quasi-equilibrium) ex¬ 
pansion or compression. In the top panel of Fig. 4, we 
plot CSDs for simulations of the u(r) potential corre¬ 
sponding to ntgt = 32 at 77 = 0.06 at various terminal 
(equilibrated via long runs between compressions) pack¬ 
ing fractions 77 = [ 0 . 02 , 0 . 12 ], where it is apparent that 
the CSDs possess primary peaks (30 < n* < 35) centered 
near the original targeted value. 

Consistent with the notion that free monomer parti¬ 
cles represent the “vapor” in a microscopic liquid-vapor 
coexistence, as 77 decreases, the integrated amount of 
monomer and other small aggregates increases (as ex¬ 
pected from energy-entropy compensation arguments) 
and the main peak in the CSD shifts left relative to the 
original position. Notably, at the largest packing frac¬ 
tions, 77 = 0.10 and 0 . 12 , secondary peaks emerge with 
local maxima at n ss 2 n*, which can be attributed to con¬ 
figurations where at least two particles from neighboring 
clusters come within r cu t of one another with some fre¬ 
quency. In fact, for 77 = 0.12, there are multiple peaks in 
the CSD at appoximate intervals of n* extending out to 
n « 6 n* clusters, where peak height is negatively corre¬ 
lated with n. 

The robustness of n* and the appearance of features 
in the CSDs at intervals of n* upon compression indicate 
that the clusters in the IC systems remain well differen¬ 
tiated and of the preferred size despite greater proxim¬ 
ity (and close contact). This is in contrast to analogous 
CSD measurements for SALR mixtures that form fluid 
(non-crystalline) clusters of similar characteristic size at 
77 = 0.06, which are shown in the bottom panel of Fig. 4. 
(The SALR parameters were chosen to produce compa¬ 
rably sized clusters, n* = 35, at 77 = 0.06). First, it 
is clear that, at the reference volume fraction 77 = 0.06, 















the SALR fluid has a considerably broader distribution 
of cluster sizes than the IC fluid, and monomer remains 
the dominant species (this is true for any given refer¬ 
ence r 7 , except when generating highly arrested percolat¬ 
ing gel states). Moreover, the n* peak is considerably 
more sensitive to changing 77 than in the IC systems, and 
for 77 > 0.06, there is an increasingly wide distribution 
of competitive cluster sizes. In other words, the clusters 
undergo continuous and unorganized growth with poor 
cluster distinguishability. 

These dichotomies in cluster size-specificity and dif¬ 
ferentiation at higher 77 between IC and SALR systems 
exhibiting amorphous clusters coincides with another be¬ 
havioral difference: the clusters in the IC systems, while 
remaining as distinguishable entities with intracluster flu¬ 
idity, self-organize at the COM level into crystalline su¬ 
perlattices for 77 > 0.08, indicating the density range 
where the liquid state of clusters becomes thermodynam¬ 
ically unfavorable (or even unstable).[ 66 ] In contrast, the 
SALR fluid remains disordered at the cluster level for all 
77 . The ability of the amorphous clusters in the IC system 
to self-assemble into a lattice strongly supports the inter¬ 
pretation of clusters as renormalized entities. Superlat¬ 
tice formation also attests to the notable monodispersity 
posessed by the clusters as the presence of polydispersity 
inhibits crystalline phases (for kinetic and/or thermody¬ 
namic reasons [67, 68 ]). 

While the 77 at which the IC clusters form superlat¬ 
tices is, at first glance, quite low, we note that the effec¬ 
tive packing fractions r] e g of the whole clusters (treating 
them as renormalized objects; see Section 2.1) are con¬ 
siderably higher. For monodisperse clusters of n* = 32, 
77 eff ~ 0.44 at 77 = 0.08 and r) e s ~ 0.55 at 77 = 0.10, condi¬ 
tions which approximately correspond to those at which 
crystallization is induced in simple fluids dominated by 
steep interparticle repulsions [67, 68 ]. Incidentally, this 
tendency toward COM crystallization makes it difficult 
to obtain convergence in the IBI scheme at similarly high 
effective packing fractions. 


3.3. Cluster persistence and particle motions 

In Fig. 5, we consider the temporal cluster persis¬ 
tence and single-particle (i.e., monomer) dynamics of IC 
(n* = 32 and SALR n* ~ 32 systems, where the FSAC 
correlation function profiles < 1 >(f) quantitatively demon¬ 
strate that the IC systems exhibit signficantly longer 
cluster “lifetimes” than their SALR counterparts. To 
wit, as shown in Fig. 5(a), the half-life values ti /2 (i.e., 
times at which $(f) = 0.5) of the IC systems are ap¬ 
proximately an order of magnitude or more greater than 
the values for SALR systems for given packing fractions. 
We consider this strong cluster fidelity a byproduct of 
generating highly monodisperse spherical clusters (and 
vice versa) like those illustrated in Fig. 2. If, instead, 
clusters are frequently exchanging constituents with each 
other or the monomer “vapor”, they will tend to exhibit 
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FIG. 4. (upper) CSDs corresponding the potential optimized 
for ntgt = 32 at 77 = 0.06 at volume fractions ranging from 
0.02 —> 0.12. (lower) CSDs calculated for an SALR potential 
(X = 5.7, Q = 0.2, £ = 2.0) that yields clusters with size 
n* fts 32 at 77 = 0.06, used for the same series of volume 
fractions. 


rather highly fluctuating (and/or instantaneously non- 
spherical) interfaces. Thus, instantaneous configurations 
will appear less superficially monodisperse and CSD pro¬ 
files will be necessarily broader. This notion is consistent 
with the CSDs in Fig. 4, where the less exchange-prone 
IC systems exhibit much greater size-specificity than the 
SALR mixtures (one can also compare the cluster snap¬ 
shots in Fig. 2 with those in Fig. 4 of a previous publi¬ 
cation [31]). 

In terms of density dependence, we observe that clus¬ 
ter persistence is negatively correlated with 77 for all n tg t 
cluster sizes considered, as typified by the nt g t = 32 re¬ 
sults shown in Fig. 5. This is easily understood by con¬ 
sidering that increasing 77 necessarily places clusters in 
close proximity, where their internal density fluctuations 
help facilitate the transfer of monomers between clus¬ 
ters. This qualitative effect is, of course, relevant for both 
IC and SALR systems, though the precise nature of the 
density-dependence differs due to the relative monodis¬ 
persity and sphericity of the IC clusters. 

However, an additional consideration when comparing 
the persistence and monodisperity of IC and SALR clus¬ 
ters is whether the SALR systems are allowed to form 
microcrystalline clusters (a phase change), as opposed 
to the amorphous clusters considered in Fig. 5(a). In 
prior work, we showed via time-lag snapshots that single- 
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FIG. 5. (a) Cluster persistence (FSAC) correlation function 
4>(t) for n * = 32 IC systems (unfilled symbols) and n* ~ 
32 SALR mixtures (filled symbols) at packing fractions 77 = 
0.02 (diamonds), 0.04 (squares), and 0.06 (triangles). The 
attractive strengths x i n the Equation 5 potentials at the 
three packing fractions are \ = 6.1, 5.9, and 5.7, respectively, 
and the repulsions are defined by Q = 0.2 and £ = 2 . 0 . (c) 
Single-particle mean-squared displacements (MSDs) for the 
same IC and SALR systems in (a). Insets (b) and (d) compare 
data for amorphous (mixture) and inicrocrystallizing (single¬ 
component) SALR models with n* « 32, where the potentials 
are defined by x = 5-4, Q = 0.2, £ = 2.0 and x = 6.0, Q = 0.5, 
£ = 2.0, respectively. Note that for visual clarity, lines trace 
all available data in each panel while symbols do not. 


component SALR models tend to form crystalline clus¬ 
ters with greater temporal persistence than slightly size- 
polydisperse (at the monomer level) mixtures that thwart 
crystallization [31]. In Fig. 5(b), this is shown quanti¬ 
tatively for n* ss 32 clusters at = 0.125, where we 
find the single-component SALR clusters exhibit fi /2 val¬ 
ues orders of magnitude longer than amorphous clusters 
of the SALR mixtures. These crystallized clusters also 
display half-lives much longer than even their IC coun¬ 


terparts: for example, recalling that increasing density 
accelerates the FSAC decay rate, the single-component 
n* « 32 SALR clusters at 77 = 0.125 exhibit ti / 2 ~ 10 4 , 
which is comparable to that of n* = 32 IC clusters at 
only 77 = 0.04. Given the orders of magnitude discrep¬ 
ancy in ti /2 between the single-component and mixture 
SALR models, it seems reasonable to ascribe the qualita¬ 
tive differences in cluster monodispersity and persistence 
to the microcrystallinity in the former, as opposed to the 
modest perturbation to the liquid state resulting from 
the weak polydispersity that we have employed. [61] 

Comparing Figs. 5(a) and 5(c), we observe that trends 
in amorphous cluster persistence are complemented by 
the single-particle mean-squared displacement (MSD) 
profiles, where particle motions comprise both intraclus¬ 
ter diffusion and slaved motion due to diffusion of en¬ 
tire clusters. First, we note that the IC systems show 
slower single-particle dynamics relative to SALR mix¬ 
tures, presumably because the slaved-motion effect per¬ 
sists out to longer timescales. Interestingly, for the IC 
systems, we also find the emergence of transient plateaus 
in the MSD profiles by 7 / = 0.06 for n tg t = 32 (and 
slightly higher 77 for smaller 7i tg t)- (A similar plateau 
also emerges for the highly persistent clusters associated 
with the single-component microcrystallizing SALR sys¬ 
tem; see Fig. 5(d).) Such a feature is typically observed 
in the context of cooperative glassy single-particle dy¬ 
namics in dense and or supercooled fluids [67, 68 ]. For 
these systems, on the other hand, it is intuitive that this 
signature is a consequence of the whole-cluster r / e g being 
much larger than rj for systems of well-defined spherical 
clusters (as described in Section 3.2). 

Thus, even for rather low packing fractions rj < 0.08, 
the presence of persistent and highly distinguishable clus¬ 
ters of an appreciable size (e.g., n* > 16) drives the emer¬ 
gence of cooperative whole-cluster dynamics, which is 
then reflected on the single-particle level. In other words, 
the 77 -range over which ICs truly diffuse around one an¬ 
other in a fluid-like manner is quite small. In contrast, 
no such shoulders in the MSDs emerge for the amor¬ 
phous SALR cluster phases up to and above 77 = 0.20 
(not shown), pointing to the quite rapid exchange and 
reformulation of clusters in non-crystallizing SALR fluids 
(provided the short-range attractions are not so strong as 
to generate dynamically arrested gel phases). 


3.4. Optimized potentials 

In Fig. 6 , we turn our attention to the pair poten¬ 
tials u(r) that result from the IBI optimization to yield 
ICs, showing particular examples for various 7it g t at 
77 = 0.04 (all 10 potentials-3 for n tg t =8 and 32, and 4 
for n tg t = 16-are provided in Supplementary Material). 
The sensitivity of the potentials to density is minimal [69] 
and, in all cases, the potential is dominated by a broad 
attractive basin terminated by a repulsive barrier that 
falls off quickly about its maximum at r rep . As demon- 
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FIG. 6. (a-c) Representative examples at g = 0.04 of the IBI 
optimized pair potentials which yield the target radial distri¬ 
bution functions for n =8, 16 and 32 respectively, (d) Op¬ 
timized potential r rep values and their corresponding target 
simulation average radii of gyration {R) values for all densi¬ 
ties studied. The linear fit (as indicated by the dashed line) 
is constrained such that r rep (0) = 0. 


strated in Fig. 6(d), the lengthscale r rep is intimately 
related to the cluster size, as it is directly proportional 
to the average cluster radius of gyration (R) extracted 
from the constrained simulations and has virtually no 
dependence on r] for a fixed cluster size. Interestingly, 
this direct proprotionality between r rep and ( R ) incor¬ 
porates the physically intuitive constraint that r rep must 
vanish when (R) = 0. This implies that the cluster sizes 
8 < n* < 32 fall in the “large-cluster” asymptotic limit 
where the discrete nature of the particles comprising the 
clusters is unimportant. 

To understand how these pair potentials strongly limit 
cluster growth to create renormalized objects, we turn 
to Fig. 7, where we show that the broad attractive wells 
drive local densification while the repulsive barriers col¬ 
lectively generate strong repulsive coronas around the ag¬ 
gregates. Fig. 7(a) shows the average potential energy 
that a particle placed in or near an n = 32 cluster expe¬ 
riences as a function of distance from the cluster COM 
using the pair potential optimized at g = 0.06. We calcu¬ 
late this interaction by averaging over particle positions 
from a simulation of an isolated cluster and either include 
or exclude the hard-core component (r < d) from the 
calculation, where the latter is done to provide a highly 
averaged measurement for the intracluster environment. 
From the former, it is clear that before any particle ap¬ 
proaches the cluster sufficiently closely to sample the very 
steep effective potential derived from excluded volume ef¬ 


fects, it “sees” an additional repulsive barrier at larger r 
that, in effect, terminates growth. From the latter, we 
see that the particles which comprise the cluster are sit¬ 
uated within a spherical attractive region by way of com- 
parision to the extent of the radial density profile of the 
cluster (Fig. 7(b)). 

In Fig. 7(c), we show a representative two-dimensional 
slice of the potential energy landscape for a single n = 32 
cluster configuration, which illustrates that the repul¬ 
sive corona is instantaneously quite strong (2 k^T < 
U(x,y) < 6kBT) and has a thickness comparable to 
the cluster radius. This latter observation can be un¬ 
derstood as follows: Fig. 7(d) shows an appropriately 
scaled schematic of a cluster that is also aligned with 
the heat map in panel (c), where we show the poten¬ 
tial due to the red-colored particle as a heat map. To 
the right, the repulsive barrier of the highlighted particle 
roughly coincides with the opposite edge of the cluster, 
building in size specificity. However, the isotropic nature 
of the pair potential necessitates that this same parti¬ 
cle also contributes repulsions in the opposite direction, 
where its repulsive barrier approximately coincides with 
the “outer” edge of the repulsive corona. This slaving 
of the outermost cluster repulsion range to cluster size 
may explain an outcome of the IBI optimizations: as 
ntgt was increased, thicker protective shells surrounding 
the clusters (as quantified by ci e ff) were required in the 
constrained simulations in order to stabilize the subse¬ 
quent IBI scheme. While not precluding the existence of 
cluster forming pair interactions that do not obey such 
size-thickness slaving; the above interpretation sugges¬ 
tive that the most “reasonable” pair interactions capable 
of generating clusters do. 

Interestingly, our pair potentials may be experimen¬ 
tally realizable via charged-monolayer protected gold 
nanoparticles, as demonstrated by Alexander-Katz and 
coworkers [70]. At the level of the pair free energy change 
between two such particles along a radial coordinate, they 
found from van der Waals, electrostatic, phobic and en¬ 
tropy contributions that a wide attractive basin followed 
by a repulsive hump can be realized. Importantly many 
tunable parameters exist in this experimentally realiz¬ 
able system making this a promising avenue for future 
research. 


3.5. Influence of long-range interactions 

While the main features of the potentials occur on the 
order of a few particle diameters in length, all of the 
optimized pair potentials also possess weak, oscillating, 
longer ranged tails. In the interest of simplifying the op¬ 
timized potentials, we consider the impact of eliminating 
these tails by cutting and shifting the ntgt = 16 opti¬ 
mized potentials at the first minimum beyond the main 
repulsive hump (as labeled in Fig. 6(c)).[71] 

In Fig. 8(a), we compare the CSDs corresponding to 
both the fully optimized and the truncated potential for 
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FIG. 7. (a) Average potential energy for interaction of a clus¬ 
ter n = 32 with a single particle with (dashed) and with¬ 
out (solid) the HC component of the IBI-optimized potential, 
(b) Radial density profile for a single cluster, (c) A two- 
dimensional slice of the potential energy corresponding to a 
representative cluster, (d) Schematic of a cluster with the 
optimized potential (n = 32) overlaid for a selected particle. 
The outer rim of the cluster (dotted circle) has a radius of 2d 
based on the intercluster radial distribution function (panel 
b). Data in panels (a-c) are generated by a simulation of a 
single cluster (n = 32) with a constrained radius of gyration, 
and the data in panels (a-b) are radially averaged. 


= 0.04: both potentials result in a clustered system 
that is fluid at both the intra- and inter-cluster level, 
indicating that the longer-ranged tail is not a strict re¬ 
quirement for IC-like behavior. However, the truncated 
potential yields clusters that are, on average, smaller 
(n* = 14) than nt g t as well as less size-specific, as ev¬ 
ident from the nearly two orders of magnitude increase 
in the CSD minumum between the n = 1 and primary 
n* peaks. Nonetheless, more ideal clustering behavior 
can be restored by decreasing temperature by less than 
0.2kBT , as shown in Fig. 8(b). These changes with T can 
be understood in the context of shifting the microscale 
liquid-gas co-existence towards the liquid side (i.e., favor¬ 
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FIG. 8. (a) and (c) Effect of cutting and shifting one ntgt = 16 
optimized potential (at the first minimum beyond the primary 
repulsion) on the CSD and RDF respectively (b) Evolution of 
the cut potential CSD from panel (a) with modest tempera¬ 
ture rescaling. 


ing clusters over monomers) upon cooling. Overall, these 
observations support the notion that while the longer- 
ranged tail modulates the cluster-to-monomer ratio (and 
thus the monodispersity of the aggregates), it is not re¬ 
quired to form ICs; thus, cluster formation in our sys¬ 
tems is predominantly a result of the competitive broad 
attraction well and repulsive hump. 

The RDFs for the full and truncated potentials are 
compared in Fig. 8(c). The obvious depression of the 
intercluster depletion region suggests that upon trunca¬ 
tion, we effectively increase the repulsive footprint of a 
cluster. As a result, intercluster crystallization will likely 
be harder to avoid. This is confirmed in Fig. 9 as the 
cut potentials associated with both of the higher density 
cases (rj = 0.06 and 0.08) formed superlaticces. 

Interestingly, the superlattice states have much bet¬ 
ter cluster definition and size preservation than the two 
lower density systems, = 0.02 and 0.04, that remained 
as IC fluids. Whether the enhanced cluster preservation 
is a byproduct of crystallization or simply higher den¬ 
sities (less void space) is not entirely clear. However, 
we do note that, in general, crystallization in systems 
of predominantly repulsive particles requires a high level 
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FIG. 9. Effect of cutting and shifting the ntgt = 16 optimized 
potentials on the CSD for all volume fractions. 

of monodispersity. This is something our emergent (not 
quenched in) cluster entities can adaptively realize in or¬ 
der to explore more thermodynamically desirable, crys¬ 
talline regions of phase space. In stark contrast, the fluid 
state does not have any natural propensity towards well- 
defined, monodisperse clusters, (polydispersity and dis¬ 
order is favored) thus making the design of an IC fluid 
all the more compelling. 


4. CONCLUSIONS 

In this paper, we demonstrated a novel application of 
standard inverse design methodology, namely, the tar¬ 
geted fabrication of liquid state structure. This approach 
was successfully applied to discover a new class of pair po¬ 
tentials that stabilize ideal cluster (IC) fluids, comprised 
of long-lived, monodisperse, spherical fluid droplets with 
good center of mass mobility. As compared to equilib¬ 
rium fluids of amorphous clusters generated via SALR 
potentials, the IC fluid states of the designed potentials 
displayed much greater size-specificity (i.e., more sharply 
peaked CSDs) with cluster sizes that were less sensitive 
to changes in overall density. Furthermore, using a new 
measure for cluster lifetime, the ICs were shown to persist 
longer than comparably sized amorphous SALR clusters 
and to maintain their identities on timescales relevant for 
cluster diffusion. 

The ability of the optimized potentials to stabilize ICs 
can be understood in terms of their general features: the 
broad (rather than short-range) attractive wells allow for 
many particles to closely pack before the relatively nar¬ 
row repulsive barriers are sampled. Moreover, the repul¬ 
sive barrier directly encodes the scale for aggregation (as 


evidenced by the proportionality of the barrier length- 
scale r rep and cluster radius ( R )) and furnishes the indi¬ 
vidual clusters with well-defined repulsive shells. By con¬ 
trast, only from a Fourier-space perspective can a pref¬ 
erential length scale for ordering be gleaned from SALR 
potentials [31]. It seems intuitive that this disparity be¬ 
tween the IBI-optimized and SALR potentials is responsi¬ 
ble for the former’s enhanced cluster size-specificity (even 
under compression). While there is a continuum of pos¬ 
sible SALR functional forms that could in principle yield 
ideal clusters, our results suggest the opposite given that 

(1) the IBI scheme did not result in SALR potentials, and 

(2) no pairwise SALR fluids have been reported that dis¬ 
play IC-like behavior for small cluster sizes (though low 
density systems of large clusters, albeit with significant 
free monomer, have been seen [35]). 

In addition to their broad attractive basins and sharp 
repulsive barriers, the optimized potentials also possessed 
weak longer-ranged oscillatory tails; however, these tails 
were found to be non-essential for IC formation. Upon 
truncating the optimized potentials beyond the repul¬ 
sive barrier, systems still displayed IC-like behavior, 
though with slightly reduced size-specificity (tending to¬ 
ward n* < ntgt)- To intentionally steer the IBI opti¬ 
mization towards shorter-ranged potentials, it may be 
fruitful to modify the constrained MC simulations to re¬ 
produce the general features of the g(r) profiles associ¬ 
ated with the truncated potentials. For instance, because 
truncating the potential enhanced the intercluster deple¬ 
tion region, one might choose (1) a stronger radius of 
gyration constraints to densify the clusters and/or (2) 
larger repulsive cluster shells to better separate the clus¬ 
ters. Optimizing for the resulting <?tgt(0 profiles may 
then naturally yield short-ranged potentials. 

In closing, our success at optimizing for IC fluids 
demonstrates the impressive flexibility of pair potentials 
for generating intricate multiscale architectures. Even 
greater flexibility could be achieved by the inverse de¬ 
sign of more complex patchy, or fully angularly depen¬ 
dent potentials; however, new schemes for solving the in¬ 
verse statistical mechanics problem must be developed. 
An interesting application of such a method would be 
to extend our current results to a patchy particle model 
where the greater degrees of freedom, while adding com¬ 
plexity, would likely allow for practically simpler (though 
likely of similar spatial range) interactions. More gener¬ 
ally, liquid-state inverse design also opens the door to 
discovery of new kinetically arrested materials: given the 
propensity for the clusters to act as renormalized objects 
(e.g., COM crystallization), this could include glasses or 
gels of supraparticles created by quenching or compres¬ 
sion out of the fluid state. 
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I. IBI POTENTIALS: ALL DENSITIES 

In Figs S1-S3 we provide the potentials optimized at each volume fraction studied in 
the main text (rj =0.02, 0.04, 0.06, and the additional case of 0.08 for n tgt = 16). Evident 
from these figures is the very minor role that density plays in determining the optimized 
potentials. Structural features largely remain fixed spatially and only a weak depression of 
the amplitude is found with increasing density. The larger amplitudes found at the lower 
densities presumably reflect the greater entropic drive the potential is competing with to 
release particles into the void region between clusters (i.e., greater space). 



FIG. SI. Optimized pair potentials for n tg t = 8. 


* truskett@che.utexas.edu 


1 








4 


a =16 

tgi 



0 2 4 6 8 

r/d 


FIG. S2. Optimized pair potentials for ntgt = 16. 



FIG. S3. Optimized pair potentials for n tg t = 32. 

II. ADDITIONAL BOND ORIENTATIONAL ORDER ANALYSIS DATA 

In Figs S4 and S5, we provide the bond-orientational order parameter distribution com¬ 
parisons for n tg t = 8 and 16, respectively. At a qualitative level, the agreement between the 
constrained and optimized systems is good; however, the agreement systematically decreases 
with decreasing cluster size as discussed in the main text. 
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FIG. S4. Bond-orientational order probability distributions for n tgt = 8 at packing fraction 
r] = 0.06. The filled grey and open blue curves indicate constrained and optimized simulations 
(using clusters with sizes between n* ± 5n as provided in the main text), respectively. 
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FIG. S5. Bond-orientational order probability distributions for ntgt = 16 at packing fraction 
r/ = 0.06. The filled grey and open green curves indicate constrained and optimized simulations 
(using clusters with sizes between n* ± 5n as provided in the main text), respectively. 
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